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ABSTRACT 

Cooling is the main process leading to the condensation of gas in the dark matter potential 
wells and consequently to star and structure formation. In a metal-free environment, the main 
available coolants are H, He, H 2 and HD; once the gas is enriched with metals, these also 
become important in defining the cooling properties of the gas. We discuss the implementation 
in Gadget-2 of molecular and metal cooling at temperatures lower that 10 4 K, following the 
time dependent properties of the gas and pollution from stellar evolution. We have checked 
the validity of our scheme comparing the results of some test runs with previous calculations 
of cosmic abundance evolution and structure formation, finding excellent agreement. We have 
also investigated the relevance of molecule and metal cooling in some specific cases, finding 
that inclusion of HD cooling results in a higher clumping factor of the gas at high redshifts, 
while metal cooling at low temperatures can have a significant impact on the formation and 
evolution of cold objects. 
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1 INTRODUCTION 

The understanding of cosmic structure formation and evolution is 
one of the most outstanding problems in astrophysics, which re- 
quires dealing with processes on very large scales, like galactic or 
cluster properties, and, at the same time, very small scales, like 
atomic behavior of gas and plasma. To join these extremes, it is 
fundamental to include atomic physics into astrophysics and cos- 
mology. In fact, only with a unified study it was and is still possible 
to justify many physical phenomena otherwise not explained, like 
for example the very well known OIII forbidden line, typical of 
many gaseo us nebu l ae. In teresting introductions into this subjec t 
are found in | Spitzerl(ll978h . the paper review by Oste rbrockl d 1988h 
and lOsterbrockUl989l) . 

Nowadays, one of the main links between "small scales" and 
"large scales" seems to be the cooling properties of the gas, as, 
to form cosmic structures, it is necessary for the gas to condense 
in the dark matter potential wells and e mit energy as radiation 
(comprehensive review s on the topic are Barkana & Loebl [20011 : 
ICiardi & Ferrarj2005T) . For this reason it is fundamental to investi- 
gate the chemical properties of molecules and atoms and their cool- 
ing capabilities. In the standard cosmological scenario for structure 
formation, the first objects are supposed to form in metal-free ha- 
los with virial temperatures lower than 10 4 K, for which atomic 
cooling is ineffective. In such physical conditions the most effi- 



cient coolants a re likely to be molecules (e.g. Lepp & Shulllll984l : 
iPuvetalJI 19931) . 

As hydrogen is the dominant element in the Universe, with a 
primordial mass fraction of about 76%, we expect that the de- 
rived molecules will play a role in the cosmologic al gas chemistry. 
The fi rst studies in th is direction were made by ISaslaw & Zipovl 
( Il967t) followed by IPeebles & Dickel Jl968h and many others 
Irlollenbach & McKeell 19791 : 1 Abel et alii 1 9971 : iGalli & Pallall998t 
IStancil et alj 19981) . who highlighted the importance of H2 in cool- 
ing gas down to temperatures of about 10 3 K. 
In addition, one should also consider that, besides hydro- 
gen, nucleosynthesis calculations predict the existence of pri- 
mordial deuterium and lithium. Rec ent measurements fr om a 
metal-poor damped Lyman a system lO' Meara et afl l2006h give 
log(D/ H) = —4.48 ± 0.06 and are consistent with o ther obser- 
vations jBurles & Tvtie3ll998l : Ipettini & Bowenll200ll) . while the 
abundance of Li (around 10 -10 ) is not very well determined and 
can vary by a factor of two or three wh en compared to the measure- 
ment s in the atmospheres of old stars dKom et alj2006l ; |Yong et all 
2006). Other molecules deriv ed from Li (e.g. LiH and LiH + ) have 
much lower abunda nces (e.g. lLepp & Shulll984l : |Puv et al.ll 1 9931: 
iGalli & Pallal 19981) . 

Another potentially interesting molecule is HD. Due to its perma- 
nent electric dipole momenQ, HD has higher rotational transition 
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Some values of the permanent HD electri c dipole moment found 
in the literature are D = 8.3 ■ 10~ 4 debye lAbgrall etaHll982l) and 
D = 8.51 ■ 10~ 4 de bye jThorson et all 19851) . The first data date back to 
iMcKellar et alj j 19761) : for a theoretical, ab initio, non relativistic, pertur- 
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probabilities and smaller rotational energy separations compared to 
H2 and thus, despit e its low abundance of log(HD/H2) ~ — 4.5 
(e.g. llxpp & Shullll 19841 : Ipuv et alJI 19931: IGa lli & Pall al 19981). HP 
can be an efficient coolan t ( Flower 200(1: IGalli & Pallal 120021 : 
iLipovka et all 120051 : lAbgrall & Rouefd l2006h and bring the gas 
in primordial halos to temperatures of the order of 10 2 K. This 
results into a smaller Jeans mass and a more efficient fragmen- 
tation process. For halos with virial temperatures in the range 
10 3 K - 10 4 K, HD cooling can be as relevant as H2, while its 
effects are expected to be minor f or larger halos (see iRipamontil 
l2007l : IShchekinov & Vasilievll2006l) . 

The formation of primordial struct ures and stars hav e been 
investigated by many authors (lik e iBromm et alj 1 19991 120021 : 
Yoshid a et alj|2006l : lKarlssodl2006h but our understanding of the 
problem is still limited, because we are lacking informations on 
all those feedback effects (such as metal pollution, mass loss and 
energy deposition from the first stars) that profoundly affect it. In 
particular, it is now commonly accepted that the presence of met- 
als, by determining the cooling (and thus fragmentation) properties 
of a gas, influences the shape of the initial mass function (IMF), 
leading to a transition from a top-heavy IMF to a Salpeter-like 
IMF, when a critical metallicity - varying b etween ~10~ 6 Z(T) and 
^10~ 3,5 Z(T), according to different authors JSchneider etai]|2003l : 



IBromm & Loebl2003l : Schneider et al.l2003) - is reached; observa 



ionally, there a r e only few constraints 1 Frebel et a l J 2007b 



ily 1 

iTornatore et al.1 d2004r) have presented the first implementation of 
a detailed chemical feedback m odel in the numerical code Gadget 
(other works on this subject are [Raited et al.ll 19961: [Gnedin 



iKawata & Gibson! 120031 : iRicotti & Ostrikers HoolhT through which 
they study metal enrichment for different feedback/IMF scenarios. 
In this study we discuss the implementation in Gadget of molecular 
and metal cooling at temperatures below 10 4 K and we present a 
scheme able to deal both with primordial and metal enriched com- 
position. In details, we extend the existi ng implementation in Gad- 
get of H2 chemistry ( Yoshi da et al. 2003), in order to include HD, 
HeH + and metal cooling at those low temperatures. Indeed, these 
species are expected to be relevant for the formation and evolution 
of cold objects. 

The paper is organized in the following way: in Section [2] we 
describe the computations of deuterium chemistry (Section l2.ll ). 
metal lines (Section 12.21 ) and their cooling capabilities (Section 
|2.3I >; in Section [3] we perform tests of our numerical implemen- 
tation about chemical abundance evolution (Section l3.lt . cosmic 
structure formation (Section 13.21 ) and cluster evolution (Section 
13. 3b ; in Section[4] we discuss the results and give our conclusions. 



2 METHODS AND TOOLS 

In the commonly adopted scenario of structure formation, objects 
form from the collapse, shock and successive condensation of gas 
into clouds having a typical mass of the order of the Jeans mass. 
This process requires the gas to cool down, i.e. the conversion of ki- 
netic energy into radiation that eventually escapes from the system. 
This can occur via inelastic collisions which induce atomic elec- 
tronic transitions to upper states, followed by de-excitations and 
emission of radiation. Details of the cooling process will depend 
on the type of elements considered and of transitions involved. 



In a standard primordial environment, the main coolants are ex- 
pected to be hydrogen, helium and some molecules like H2 and 
HD; if the medium is metal enriched, the heavier elements become 
important coolants, thanks to a larger number of possible atomic 
transitions with different energy separations. 
The relevant quantity describing the cooling properties of a plasma 
is the energy emitted per unit time and volume, i.e. the cool- 
ing function (we will indicate it with A, adopting c.g.s. units, 
erg cm -3 s _1 fl 

The characteristic time scale for the cooling, determined by A, is 
important to discriminate whether the gas can cool during the in- 
fall phase in the dark matter gravitational potential wells: structures 
are able to form only if the cooling time is short enough compared 
to the free-fall time. 

In the present work, we focus on the effects of molecules and metals 
in gas at low temp eratures. In particular , we will include their treat- 
ment in Gadget-2 jSpringel et al. 2001; Springe] 2005). This code 
uses a tree-particle-mesh algorithm to compute the gravitational 
forces and implements a smoothed-particle-hydrodynamics (SPH) 
algorithm to treat the baryons. Moreover, it is possible to follow the 
main non-equilibrium reactions involving electrons, h ydrogen, he- 
lium and H2, with the respective ionization states dYoshida et al.l 
120031) . Stellar feedback processes and metal release from SNII, 
SNIa and AGB stars are also included together with metal cool- 
ing at temperatures h igher than 10 4 K (for a detailed discussion see 
ITornatore et al.ll20071 . and references therein). 
In the following, we are going to discuss in detail our HD and metal 
line treatment at T < 10 4 K. 



2.1 HD treatment 

The HD molecule primarily forms through reactions between pri- 
mordial deuterium and hydrogen atoms or molecules: a com- 
plete model for the evolutio n of HD involves 18 reactions 
dNakamura & Umemural |2002|) . but, as their solution becomes 
quite computationally expensive when implemented in cosmolog- 
ical simulations, w e use only the set of reactions selected by 
IGalli & Pall j §002), which are the most relevant for HD evolu- 
tion: 



H 2 + D -> HD + H 
H 2 + D + -> HD + H+ 

which lead to HD formation; 

HD + H -> D + H 2 
HD + H+ -» D+ + H 2 

for HD dissociation and H2 formation; and 

H++D -» H + D+ 
H + D+ -» H + + D 



(1) 

(2) 

(3) 
(4) 

(5) 
(6) 



for charge exchange reactions. 

From reactions Q} - we see that HD abundance primarily de- 
pends on the amount of primordial deuterium and on the H2 frac- 
tion. 

For each species i, the variation in time of its number density 

m is 



bative treatment, via radial Schroedinger equation, see also lFord & Browne] 
\ 19771) and references therein. 



Sometimes it is possible to find the same notation A for the cooling rate 
in erg cm 3 s — 1 . 



Metals and molecules in structure formation 3 



Table 1. Set of reactions in the code. 




Figure 1. Temperature evolution of the reaction rates for deuterium chem- 
istry. The labels refer to the number of the equation in the text. 
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where fc P9l i in the first term on the right-hand side is the creation 
rate from species p and q, and ku is the destruction rate from in- 
teractions of the species i with the species I; they are temperature 
dependent and are usually expressed in [cm 3 s -1 ]. 
A plot of the rates as a function of the temperature is given in Figure 
[T] and the exact expressions and references in Appendix [A] From 
the figure, it is clear that the most important reactions in the relevant 
range of temperatures are (O and {6]l, and that the HD creation rates 
of reactions (TJ and ([2} are always higher than the corresponding 
destruction rates of reactions ([3} and ©, respectively. 

We have also considered HeH + molecule evolution and 
found negligible effects on the cooling properties of the gas. The 
rates for HeH + formation and destruction are given in AppendixlAl 



We impl e ment our chemistry model extending t he code by 
lYoshida et al.l | |2003|) . which adopts the rates from lAbel et all 
( 1997), and modify it for self-consistency to obtain a set of reac- 
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P7 1 = IPeterson et all |l97ll); K AH79 = iKarpas et alj Jl979t): 
' _1982h : SK87 = IShapiro & Kai5 
GP98 = iGalli & Pallal Jl998l 



= Roberse & Dalsarno 
A97 = 



lAbel et al.l 
Standi et al.l 



ST99 



; W S02 = IWane & Standi J2002h ; 

= iGlover & Brand (2003); S04 
Yoshidaetaljfe006l) . 



S02 



Di = >J kuTii. 



Stibbe & Tennyson 
= ISavinl |2002 j 
ISavin et alj 12004 



(10) 



The number density, n\, is then updated from equation (8): 

t+At Cl +At At + nl 
Ul 1 + D l + At At ' 

We apply this treatment to all chemical species. 



(11) 



HD, HeH + (a complete list of the reactions is given in TableQJ. 
The set of differential equations is evaluated via sim- 
ple linearization, ac cording to a backward difference formula 
(Anni nos et "all ll997b : given the time step At, at each time t and 
for each species i, equation Q can be re-written as 



t+At 



At 



C: 



t+At 



D 



t+At t+At 
71 4 



(8) 



2.2 Metal treatment at T < 10 4 K 

For our calculations, we consider oxygen, carbon, silicon and iron, 
because they are the most abundant heavy atoms released during 
stellar evolution and, therefore, they play the most important role 
in chemical enrichment and cooling: indeed, supernovae type II 
(SNII) expel mostly oxy gen and carbon, while supernovae type la 
(SNIa) silicon and iron iThielemann et al.ll200ll : IPark et alJ^OoS: 



where we have introduced the creation coefficient for the species i, iBorkowski et alj2004 iMevnet et alj200fj) . 



m [cm 3 s 1 ], as 



and the destruction coefficient, in [s 1 ], as 



We make the common assumption that carbon, silicon and iron 
are completely ionized, while oxygen is neutral. This is justified 
(9) because, in a cosmological context, UV radiation below 13.6 eV 
(from various astrophysical sources, like quasars, stars, etc.) 
can escape absorption by neutral hydrogen and generate a UV 
background that can ionize atoms with first ionization potential 
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lower than 13.6 eV (like carbon, silicon and iron), while oxygen 
remains predominantly neu tral since its first ionization potential of 
13.62 eV is higher (see also lBromm & Loebl20 03: Sant oro & Sh"ul3 
l2006h . 

As in the low density regime of interest here thermodynamic 
equilibrium is never reached (see discussion of eq.U6b. the Boltz- 
mann distribution for the population of atomic levels can not be 
used. Thus, we will use the detailed balancing principle instead. 
For each level i of a given species, we impose that the number of 
transitions to that level (which populate it), per unit time and vol- 
ume equals the number of transitions from the same level i to other 
levels (which de — populate it), per unit time and volume: 



ni Pa = S Pn j P ji (ij^j). 



3 



j 



(12) 



In formula J 1 2b . Pij is the probability per unit time of the transition 
i — > j and rii and rij are the number densities of atoms in the i-th 
and j-th (with i 7^ j) level. The left-hand side of the previous equa- 
tion refers to de-populations of the i-fh level, while the right-hand 
side refers to the transitions which can populate it. 
The probability of a given transition can be easily computed once 
the Einstein coefficients and the collisional rates are known. 
The further constraint which must be satisfied is the number parti- 
cle conservation: 



rii = titot 



(13) 



j 

where ntot is the total number density of the species considered 
and rij the population of the generic level j. 

In case of collisional events, the rate at which the transition 
i — * j occurs is by definition: 



nin x jij = n,in x (uoij) — nin x / uoijf(u)d u 



(14) 



where atj is the cross section of the process, f(u)d 3 u is the veloc- 
ity distribution function of the particles (typically a Maxwellian), 
7ij is the collisional rate, rii the number density of the particles in 
the i-lh level and n x is the colliding particle number density. 
The relation between 7^ and 7ji is: 

9Hij = 9j7]ie~ PAEj \ (15) 

where gi and g 3 are the level multiplicities, /3 = (feT)" 1 , AEji 
is the energy level separation and i < j. 

In addition to collisionally induced transitions, spontaneous transi- 
tions can take place with an emission rate given by the Einstein A 
coefficient. 

It is convenient to define the critical number density for the transi- 
tion i — > j as 

n cr ,ij = — ■ (16) 
Hi 



This determines the minimum density above which thermal 
equilibrium can be assumed and low density deviations from the 
Boltzmann distribution become irrelevant. At densities below 
n cr ,ij, we expect values of the excited level populations lower 
than in the thermodynamic limit, because of the reduced number 
of interaction^. 



For a two-level system, the low density level populations aris- 
ing from electron and hydrogen impact excitations can be found by 
solving the system of equations resulting from conditions U3\ and 
G3: 



f Til + Tl2 = Titot 



tninH7!|+nin e 7i2 -n 2 n H 7|f -n 2 n e 7Si -ri2^2i =0 



(17) 



where uh and n e are the hydrogen and electron number density, 
while 7^ and 712 are the H-impact and e-impact excitation rate. 
The solution of l !17b is: 



ni _ 72i + l2in e /n H + A 2 i/n H 

ntot 75 + J21 + (7i2 + 7fi) n e /n H + A 2 i/n H 

112 _ Tig + li2ne/n H 

ntot 7(2 + I21 + (7i2 + Tli) n e /n H + A 2 i/n H ' 
The ratio between the two level populations 



7i2 +lt 2 n e /n H 



n2 

ni 721 + 7f 1 n <= / nH + A21/TIH 



712 



72i + A 2 \jnH 



(18) 
(19) 

(20) 
(21) 



will in general deviate from the Boltzmann statistic, because 
the spontaneous emission term dominates over the collisional 
term at low densities. In a neutral dense gas, instead, the 
level population saturates and simply reduces to a Boltzmann dis- 
tribution, independently from the colliding particle number density. 

In case of n— level systems, one must solve the n x n popu- 
lation matrix consisting of n — 1 independent balancing equations 
J 1 2b and the constraint of particle conservation J 1 3b . 
In the modeling, we approximate CII and Sill as a two-level sys - 
tem, and 01 and Fell as a five-level system JSantoro & S hull 2006). 
Further details on the atomic data and structures are reported in Ap- 
pendix [B] 



2.3 Cooling 

In addition to calculating the chemical evolution of the gas, we need 
to evaluate the cooling induced by different species. In the original 
code, hydrogen and helium cooling from collision al ionization, ex- 
citation and recombination i tHui & Gnedinl 19971) , Compton cool- 
ing/heating and Bremsstrahlung 1 Black 198l[) are evaluat e d. For 
the H2 and cooling, the rates quoted in lGalli & Pallal i 19981) 
are ad opted. We take the HD cooling function from lLipovka et al.l 
(2005), who consider the HD ro- vibrational structure and perform 
calculations for J ^ 8 rotational levels and v — 0, 1, 2, 3 vibra- 
tional levels. Th eir results are somehow more ac curate than other 
approximations (Flower 2000 ; G alli & Pallall 20021) and valid for a 
wide range of number densities (up to 10 4 cm~ 3 ) and temperatures 
(10 2 K - 2 ■ 10 4 K). 

In Figure|2] we show cooling functions for H2, HD, H^ molecules; 
for the latter case we distinguish between neutral hydrogen impact 
and electron impact cooling; we have assumed fractions xhd = 
10" 8 ,xh 2 = 10~ 5 ,x H+ = 10" 13 ,x c - = 10" 4 and a total hydro- 
gen number density nH = 1 cm -3 . Due to its very low abundance, 



3 The critical number density depends on the particular line transition con- 



sidered; typical values for the fine structure transitions we are mostly inter- 
ested in are of the order ~ 10 5 cm — 3 . 
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Figure 2. Cooling functions for a primordial gas with a hydrogen number 
density of 1 cm -3 and the following fractions for the different species: 

x HD = icr 8 , x H , = icr 5 , x„+ = i0" 13 , x c _ = icr 4 . The h 2 

2 

cooling function (long-dashed line) is plotted together with the HD (solid), 
H-impact (dotted line) and e-impact (short-dashed line) cooling 
functions. The fits for cooling are appropriate only for T < 10 4 K. 



Figure 3. Cooling due to metals as a function of temperature. The compu- 
tations are done for a gas with total number density of 1 cm -3 ; for each 
metal species we assume a number density of 10 -6 cm -3 and we set the 
free electron over hydrogen fraction to a value of 10 — 4 . 



Hj is less effective than neutral H2 and HD, which remain the only 
relevant coolants over the plotted range of temperature. 
We highlight that the contribution of HD to gas cooling at low tem- 
peratures is dominant in the case considered here, but its relevance 
strongly depends on the relative abundances of the species. 

The cooling for metal line transitions is computed as follows. 
In case of two-level systems, we define 



A = n 2 A 21 AE 2 



(22) 



where n 2 is the atomic excited state number density, A21 is the 
probability per unit time of the transition 2 — > 1 and AE 2 i is the 
energy separation of the levels. 

Combining ( 122b and d!9b one can write the previous equation 
as a function only of the total number density of the species 



A= 



7l2 +7l2"e/nH 



7S+72i+(7i2-h'Ii) n e /n H +A 21 /n H 



n to tA 2 iAE 21 . 



(23) 



For n e < C tlh, the previous for mula is consistent with the one 
quoted in ISantoro & Shulll §006), who do not consider electron 
impact excitation effects. Using equations d20b and l !16b , A can also 
be written as a function of the fundamental level population 



nin H Ji 2 + nire e 7i2 
n a/ng 21 + n e /n e cr21 + 1 



AE-n 



n H/nf r:21 + 1 



AE 21 



(24) 
(25) 



being n cr 21 and n^ r 21 the critical density for the transition 2 — > 1 
due to H- and e-impact excitations. 

In particular, in the low density limit (nn,e <C n cr ), the above 
equation becomes 



A ~ [n-in H ^? 2 + nin e 7f 2 ] AE 2 i 
~ nwujisAEai. 



(26) 
(27) 



In this regime, each excitation - see formulae d26b and d27b - is 
statistically followed by emission of radiation - see the general def- 
inition 

122). 

In the high density limit, one finds the expected thermodynamic 
equilibrium cooling rate 



A21 



1 + 

A21 



AE 21 + 



<?2 -3AE 21 

9i 



ile-^m — e AE 21 

1 + ^721/^721 



9l 

H g2 c ~pAE 2 i 
91 



niA 2 iAE 21 . 



(28) 
(29) 



In the right-hand side, it is easy to recognize the Boltzmann 
distribution of populations for n 2 . It is interesting to note that the 
cooling function does not depend any more on the number density 
of the colliding particles, but only on the species abundance, 
in contrast with the low density regime, where there is a linear 
dependence on both densities. 

These arguments ensure that it is safe to use formula d23b to 
compute the gas cooling for two-level atoms. 

For n— level systems, the cooling function is simply the sum 
of all the contributions from each transition 



A i 



ruAijAEi. 

i^l 0^j<i 



(30) 



In general, once the number density of the cooling species is fixed, 
we expect the cooling function to grow linearly with the colliding 
particle number density and eventually to saturate, converging to 
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Figure 4. Total cooling due to hydrogen, helium, metals, H2 and HD 
molecules as function of temperature, for gas having a hydrogen number 
density of 1 cm -3 . The fraction of H2 and HD are fixed to 10 — 5 and 10 — 8 , 
respectively. The labels in the plot refer to different amount of metals, for 
individual metal number fractions of 10 — 3 (solid line), 10 — 4 (long-dashed 
line), 10 -5 (short-dashed line) and 10 -6 (dotted line). 



the Boltzmann statistic, when the critical densities are reached. We 
see that CII, Sill, Fell saturate when the colliding particle number 
density achieves values around 10 4 cm" 3 — 10 5 cm~ 3 , while 
for OI we will have a double phase of saturation: the first one at 
~ 10 5 cm -3 involving the lower three states and the second one at 
~ 10 11 cm" 3 involving the higher two states. 

As an example, in Figure [3] we show the cooling functions 
for a total number density 1 cm -3 and for each metal species 
lCP 6 cm~ 3 ; the ratio between free electrons and hydrogen is 
chosen to be 10 -4 . With these values, the presence of electrons 
can affect the results up to 10% with respect to the zero electron 
fraction case. We also notice that all the metals contribute with 
similar importance to the total cooling function and the main 
difference in the cooling properties of the gas will depend on their 
detailed chemical composition. 

We also plot the cooling functions for all the temperature 
regime we are interested in: at temperatures hig her than 10 4 K, we 
interp olate the Sutherland and Dopita tables ( Sutherland & Dopita 
1993), at lower temperatures, we include metals and molecules 
as discussed previously. Figure [4] shows the cooling function for 
different individual metal number fractions with abundances in 
the range 1CT 6 - 10~ 3 and H 2 and HD fractions of 10~ 5 and 
10 -8 , respectively. These values for H2 and HD are fairly typi- 
ca l for the IGM gas a t the mean density (see also the conclusions 
of iGalli & Pallalll998t and references therein). In the temperature 
range 10 4 K — 10 5 K, the double peak due to hydrogen and he- 
lium collisional excitations is evident at low metallicity, while it is 
washed out by the contribution of different metal ionization pro- 
cesses as the metallicity increases. For example, complete colli- 
sional ionization of carbon and oxygen produces the twin peak at 
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Figure 5. Abundances as a function of redshift. The solid lines refer to the 
abundance evolution in a flat cold dark matter universe with h = 0.67, 
f2om = 1, ^ub = 0.037; the dotted lines refers to H2 and HD evolution 
in a ACDM model with h = 0.73, C m = 0.237, U 0A = 0.763, 
tt 0b = 0.041. 



10 5 K, while complete ionization of iron is evident at about 10 7 K. 
At temperatures lower than 10 4 K and metal fractions lower than 
~ 10 -6 , the dominant cooling is given by molecules; instead, for 
larger metal fractions the effects of metals became dominant. 
The general conclusion is that at very high redshift, when metals 
are not present, only H2 and HD can be useful to cool the gas down 
to some 10 2 K, while after the first stars explode, ejecting heavy 
elements into the surrounding medium, metals quickly become the 
most efficient coolants. 



3 TESTS 

In this section, we are going to test the implementation of HD and 
metal cooling using different kind of simulations. In particular, we 
focus on the analysis of abundance redshift evolution, cosmic struc- 
ture formation and clusters. 



3.1 Abundance redshift evolution 

As a first test, we investigate the behavior of a plasma of primor- 
dial chemical composition (i.e. with no metals) looking at the red- 
shift evolution o f the single abundanc es. Our goal is to reproduce 
the results from IGalli & Pallal fl998), who calculate the redshift 
evolution of a metal-free gas at the mean density by following 
a detailed chemical network. For this reason, here, we perform 
our non-equilibrium computations on isolated particles, includ- 
ing the following chemical species: e~, H, H + , He, He + , He ++ , 
H2, Hj, H , D, D + , HD, HeH + and assuming a flat cosmology 
with no dark energy content (matter density parameter fiom = 1), 
baryon density parameter fiob = 0.037, Hubble constant, in units 
of 100km s -1 Mpc -1 , h = 0.67 and initial gas temperature of 
1000 K. 

The evolution of the number fractions for the different 
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species is plotted in Figure |5j the electron abundance is given 
from charge conservation of neutral plasma and is normally very 
close to the H + value, this being the do minant ion. These re - 
sults are in very good agreement with those of lGalli & PallaUl998l) . 

In our set of reactions, due to the low initial gas tempera- 
ture, the collisions are inefficient to ionize helium. The inclusion 
of HeH + creation 



He + H + -> HeH + 7 

contributes to rise H 4 / abundance mainly via reaction 
HeH 4 " + H -> He + H 4 / 

and weakly decrease the H~ number fraction via 
HeH 4 " + 7 -» He + H + 



H T + H" 



2H 



(31) 



(32) 



(33) 
(34) 



where 7 indicates the photons. Because of the very low HeH 4 " 
abundance reached, there is no substantial He atom abundance evo- 
lution. 

Another caveat to take into account is the lack of reactions 
between D + and free electrons which would destroy the deuterium 
ions more efficiently, but without altering significantly the global 
amount of HD formed. We notice also the exponential decay of 
D + due to the rate coefficient of equation §5$ and the freezing out 
of H + , H2, D and HD number fractions. 

As a comparison, we also plot (dotted lines) the H2 and HD 
abundance evolution in a flat ACDM model having h = 0.73, 
fiom = 0.237, n A = 0.763, n oh = 0.041 (Spergel et al. 2006). 
The slight increment observed is due to the fact that in the cold 
dark matter cosmology the baryon fraction is about 4%, making 
the interactions among different species rarer than in the ACDM 
model, for which the baryon fraction is about 17%. In addition, 
the cosmological constant is dominant only at redshifts below one. 
The evolution of the other species is similar in both cosmologies. 



3.2 Cosmic structure formation 

To test the behavior of the code in simulations of structure 
formation and evolution and the impact of HD, we run a cosmo- 
logical simulation with the same properties and cosmology as in 
iMaio et alj ( 120061) . The main difference here is the addition of HD 
chemistry. We adopt the concordance ACDM model with h = 0.7, 
fiom = 0.3, f2ob = 0.04, floA = 0.7; the power spectrum is 
normalized assuming a mass variance in a 8 Mpc/h radius sphere 
erg = 0.9 and the spectral index is chosen to be n = 1. We sample 
the cosmological field (in a periodic box of 1 Mpc comoving 
side length) with 324 3 dark matter particles and the same number 
of gas particles, having a mass of about 1040 Mq and 160 Mq, 
respectively. The comoving Plummer-equivalent gravitational 
softening length is fixed to 0.143 kpc. This allows to resolve halos 
with mass of about 10 5 Mq. The simulation starts at z — 90 and is 
stopped at z ~21. 

We include the reactions involving e~, H, H 4 , He, He 4 , He 44 , 
H 2 , Hj, H", D, D + , HD (here we neglect HeH 4 , as it has not 
signific ant effects on the s imulation) and compare the results with 
those of lMaio etal] ( l2006l) . whose ACDM simulation has the same 
features, but the chemical set does not follow the evolution of D, 
D 4 " and HD and does not include H \ cooling. 
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Figure 6. Gas clumping factor as a function of redshift for two ACDM mod- 
els with different chemical composition. The squares refer to the clumping 
factor computed with standard atomic line cooling and H2 cooling, while 
the triangles refer to a case which includes also HD cooling. The shaded 
regions conespond to the variation of the maximum overdensity between 
100 (lower line in both cases) and 500 (upper line in both cases). 



To quantify the differences between the two runs and the effi- 
ciency of the HD cooling we calculate the gas clumping factor, C, 
in the simulation box, in the following way 



G = 



(35) 



where for each SPH particle, i, we indicate with m.i its mass and 
with pi its mass density; the indices run over all the gas particles. 
For the sake of comparison, we calculate C using only particles 
with density below a given overdensity threshold, 5m > and we make 
5m vary in the range [100, 500]. 

The results are plotted in Figure|6]for both simulations. We see that 
the inclusion of HD makes the clumping factor increase at all red- 
shifts, almost independently from the density threshold. This means 
that the gas is, on average, denser and more clumped, with an in- 
crement of about 10% at redshift ~22. 



3.3 Cluster 

So far, we have assumed either primordial gas with no metal pollu- 
tion (previous test case) or a pre-defined metallicity to demonstrate 
the effect of the presence of metals on the cooling function at low 
temperature, as in Section |2~31 Now, we are going to couple our 
cooling function with a model for the chemical enrichment and test 
this implementation within a simulation that follows the formation 
of a cluster. In addition to testing the validity of our implementa- 
tion, although there are no significant changes for the intra-cluster 
medium (ICM) to be expected, it is of interest to check whether 
there are regions inside the simulations where the polluted medium 

is cooling below 10 4 K due to its metal content 

The "zoomed initial condition technique" dTormen et alj 1 19971) 
is used to extract from a dark matter-only simulation with 
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Figure 7. Distribution of the particles of a cluster simulation in the T-A 
space. The hot and thin intra-cluster medium populates the bottom right 
area. The particles within collapsed objects, which represent very dense re- 
gions of the simulation, populate the high temperature cooling function (up- 
per branch). In addition, metal enriched particles undergo metal line cooling 
in the low temperature regime (below ~ 10 4 K). The three-pointed star 
symbols correspond to particles which are located within twice the virial 
radius of the cluster and have a temperature lower than 8000 K. 




10= 10= to 31 10= 10= 10= 10= 10= 10= 10= 10= 
p [g cm" 3 ] 

Figure 8. Distribution of the particles of a cluster simulation in the phase 
diagram. The hot and thin intra-cluster medium populates the central-left 
area of the plot, while the dense and cooled regions of the simulation are 
represented in the lower-right part. Particles heated by feedback effects 
are represented by points in the central-right side. The three-pointed star 
symbols correspond to metal enriched particles which are located within 
twice the virial radius of the cluster and have a temperature lower than 
8000 K. We remind that the critical density of the Universe at the present 
is p ,cr =! 1.9 ■ 10~ 29 fc 2 gem,- 3 . 



box size of 479Mpc/h (we adopt a ACDM cosmology with 
H = 70 km s" 1 Mpc -1 , eg = 0.9, Qoa = 0.7, fi 0m = 0.3, 
f^ob = 0.04) a smaller region and to re-simulated it at higher 
resolution introducing also gas particles. The cluster evolution is 
simulated with about 2 ■ 10 5 particles. The comoving Plummer- 
equivalent gravitational softening length is 5kpc/h. At redshift 
zero, the selected cluster has a virial mass of about 10 14 M Q /h, 
a virial radius of about IMpc/h and a vir ial temperature of 2 • 10 7 
K (for more details see lDolag et alj|20040 . 

We start the simulation with no metallicity c ontent. Then, the meta l 
abundances are consistently derived (as in iTornatore et al.ll2007t) 
following the star formation history o f the system, accounting for 
the lifetime of stars of different mass jPadovani & Matteuccil 19931) 
distributed according to a Salpeter I MF and adopting appropri- 
ated stellar yields: we use those from Wooslev & Weaver! J 19951) 
for massive stars (SNII). lvan den Hoek & Gr oenewegen dl997l) for 
low- and intermediate-mass stars and lThielemann et alj d2003l) for 
SNIa. The underlying sub-resolu tion model for star formatio n in 
multi-phase interstellar medium ( Springel & Hernquist 2003) in- 
cludes a phenomenological model for feedback from galactic ejecta 
powered by the SNII explosions, where we have chosen the wind 
velocity to be 480km s" 1 . As we are only interested to test the ef- 
fect of the metals, we exclude H2, HD and HeH + chemistry and 
consider only atomic cooling from collisional excitations of hydro- 
gen and helium. Once the medium gets polluted with metals, their 
contribution is added. For the metal cooling of the gas abo ve 10 4 K, 
Sutherland and Dopita tables ( [Sutherland & Dopital 19931) are used. 
At lower temperatures, the fine structure transitions from OI, CH, 
Sill, Fell are included as discussed in the previous Sections. 

Figure [7J shows the cooling diagram of our simulation at red- 
shift z — 0; each SPH particle is represented by a point. In the 



plot, different areas can be identified. The one at high temperatures 
(bottom right) represents the hot ICM. When the ICM starts to get 
denser, cooling gets more efficient: the corresponding gas particles 
are represented by the points belonging to the upper branch of the 
cooling function and they are brought to lower and lower temper- 
atures. Feedback from the star formation partially pushes some of 
them away from the cooling curve to slightly higher temperatures. 
Below 10 K, only particles which are metal enriched can further 
cool down to about 300 K, while gas particles with primordial com- 
position are stacked at 10 4 K. The three-pointed stars refer to par- 
ticles within twice the virial radius and with a temperature below 
8000 K, indicating that this region of the T-A space is also popu- 
lated by gas associated with the galaxy cluster. 
The corresponding phase diagram is shown in Figure [8] The three- 
pointed star symbols are the same as in Figure [JJ The hot and 
thin intra-cluster medium populates the central-left area of the plot, 
while the dense and cool regions occupy the lower-right part. Par- 
ticles heated by feedback are represented by points in the central- 
right side (p > 10" 24 gem -3 , T > 10 5 K ). The main effect 
of our metal cooling implementation is to lower the temperature of 
the dense medium, generating the sharp triangular area visible in 
the p -T space, at T < 10 4 K and p > 10~ 26 g cm" 3 . The points 
at very low densities are associated with diffuse metal free gas; 
this suggests that the spread in the cooling diagram of Figure [JJ at 
temperature lower than 10 4 K is mainly due to different fractional 
metal enrichment of the particles, rather then to their different den- 
sities. 

As already mentioned, global properties of the ICM and star for- 
mation are not significantly changed compared to the reference 
run without the metal line cooling from fine structure transitions. 
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This happens because the simulation was merely meant to be a test 
case of the implementation of metal line cooling below 10 4 K un- 
der realistic conditions, but the halos resolved are large enough to 
cool and form stars without the aid of such cooling. In order to 
investigate in more detail the effects of the additional cooling by 
molecules and metals at low temperatures on the ICM and the star 
formation, higher resolution simulations are needed. However, this 
opens interesting grounds for further investigations on the interplay 
between formation of small objects, with virial temperatures in the 
range of interest for our extended cooling function, and metal pol- 
lution from first stars. 



4 CONCLUSIONS 

In order to understand structure formation and evolution, a detailed 
study of the chemical and cooling properties of baryonic matter is 
needed. 

In this paper, we have presented time dependent calculations of the 
cooling properties of a gas in a "low temperature" regime, using 
the contributions of several chemical species and we have tested 
the effect on cosmic structure evolution. 

Hydrogen derived molecules are effective in cooling metal-free 
gas below a temperature of ~ 10 4 K, the typical temperature range 
of primordial objects. On the other hand, when the medium is 
polluted by material expelled from stars (via SN explosions, mass 
losses in AGB phase and winds), metals are expected to become 
the main coolants. 

For these rea sons, we have extend ed previous "non-equilibrium" 
calculations dYoshida et al.l 120031) in order to include in th e 
numerical code Gadget-2 dSpringel etafl 1200 it ISpringell l2005h . 
the deuterium chemistry and follow the formation/destruction of 
HD molecule. This, together with molecular hydrogen, is able 
to cool down the gas at T < 10 4 K. Thanks to its permanent 
electric dipole mome nt, HD could allow cooling even below 10 2 K 
dYoshida et alj2006l) . 

Other molecules are not very significant for the gas cooling 
properties. 

The treatment of metal cooling at T ^ 10 4 K is inc luded using 
the tables provided by Sutherland & Dopital dl993l) . while the 
contribution from fine structure transitions of oxygen, carbon, 
silicon and iron at T < 10 4 K has been included by computing 
the populations of the levels, for each species, using the detailed 
balancing principle. More in particular, we have assumed that 
the UV radiation coming from the parent star ionizes carbon, 
silicon and iron, while oxygen remains neutral as its first ionization 
potential is higher than 13.6 eV. We deal with the gas radiative 
losses computing the detailed balancing populations of the levels 
due to collisional excitations arising from hydrogen and electron 
impacts. The cooling follows the level de-excitations. The electron 
impact excitations are also included, as a residual electron fraction 
of about 10 -4 survives in the post-recombination epoch and higher 
values are reached during the reionization process. 
On the whole, we are now able to follow the evolution of e~, H, 
H + , He, He+, He ++ , H 2 , H+, H~, D, D+, HD, HeH+, O, C+, 
Si + , Fe + ; so, the code is suitable to deal both with primordial and 
metal enriched gas. 

We have checked the validity of our scheme by comparing the 
results of some test runs with previous calculations of molecule 
abundance evolution, finding excellent agreement. 
We have also investigated the relevance of HD and metal cooling 
in some specific cases. 



Adding the deuterium chemistry and HD contribution to the 
cooling function in simulations of structure formation results in a 
higher clumping factor of the gas, i.e. clouds are slightly denser 
and more compact, at high redshifts, with respect to the case when 
only H, He and H2 cooling is considered. The difference is about 
10% at z ~ 22. 

For what concerns the role of metal cooling at T < 10 4 K, we 
have shown that their presence is relevant in this temperature 
regime. In particular, in the cluster simulations we have run, fine 
structure transitions can actually cool the local temperature down 
to 10 2 K - 10 3 K. 

In conclusion, we have implemented in Gadget-2, the most 
relevant features of gas cooling, in both pristine and polluted envi- 
ronments, for the temperature range 2.7 K — 10 9 K. We find that 
HD cooling has some influence on the high redshift gas clumping 
properties, while low temperature metal cooling has a significant 
impact on the formation and evolution of cold objects. In addition 
to investigating the above topics, this implementation can be used 
to study the detailed enrichment history of the IGM and its possi- 
ble interplay with the transition between a primordial, massive star 
formation mode and a more standard one. 
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APPENDIX A: CHEMICAL RATES 

We consider the following set of equations involving HD creation 



and destruction: 

H 2 + D -> HD + H (Al) 

H 2 + D + -> HD + H + (A2) 

HD + H -> D + H 2 (A3) 

HD + H+ -f D + + H 2 (A4) 

H+ + D -f H + D+ (A5) 

H + D + -f H+ + D. (A6) 

We use the rate coefficients from Wang & St ancill d2002l) : 

fc HD ,i = 9.0 ■ 10" 11 e" 3876/T cm 3 s _1 (A7) 

fe D ,2 = 1.6 ■ 10" 9 cm 3 s _1 (A8) 

IStancil etalJdl998h : 

fc HD ,3 = 3.2.10- 11 e- 3624/T cm 3 s- 1 (A9) 

k WA = lO^e-^cmV 1 (A10) 
and lSavinld200l) : 



7 010-10^)0.402 -37.1/T 0110-17^1.48 3 -1 /a ,in 

khd,5 = 2-10 1 e ' —3.31-10 1 cm s (All) 
fcHD,6 = 2.0610" 10 T 0396 e" 33 0/:z +-2.0310" 9 r" ' 332 cm 3 s""( L A12) 

We consider the main equations for HeH + formation and evolution 
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jGalli & Pallall 19981) . namely: 

He + H+ -» HeH+ + 7 
HeH + + H -» He + 
HeH+ +7 -» He + H+ 

and the rates from lRoberge & Dalgarnol dl982h : 

r 7.6- 10 _18 T-°- 5 cm 3 s- 1 , T < 10 3 K 

^HcH+.l ~ \ 

{ 3.45- 10" 16 T" lo6 cm 3 s- 1 , T>10 3 K 
lKarpasetal1 ( ll979h : 
k HeH + 2 = 9-1 • 10 -10 cm 3 s _1 
and lRoberge & Dalgarnol dl982h : 

fc HcH +,3 = 6.8 • 10 _1 T r 15 e 



lrpl.5^-22750/T r -1 



(A13) 
(A14) 
(A15) 



(A16) 



(A17) 



(A 18) 



In the previous expressions, T stands for the gas temperature and 
T r for the radiation temperature. 



APPENDIX B: ATOMIC DATA 

In the following, the atomic data a dopted in this paper are pro- 
vided (see also lOsterbrockl 1 19891 ; IHollenbach & McKed 1 19891 ; 
ISantoro & Shullll2006l) . We will use the usual spectroscopic nota- 
tion for many electron atoms: S is the total electronic spin quantum 
operator, L the total electronic orbital angular momentum operator 
and J = L + S the sum operator; S, L, J, are the respective quantum 
numbers and X indicates the orbitals S, P, D, F, according 
to L=0, 1, 2, 3, respectively; then 2S+1 Xj will indicate the 
atomic orbital X, with spin quantum number S and total angular 
momentum quantum number J; its multiplicity is equal to 2 J + 1. 
In the following, we are going to discuss the models adopted for 
each species and the lines considered in a more detailed way. 
We will often use the notation T100 = T/100 K. 



* CII: we model CII as a two-level system considering the fine 
structure transition (2p)[ 2 P 3/ / 2 — 2 Pi/2] between the quantum 
number J = 3/2 and J = 1/2 states. Data were taken from 
IHollenbach & McKed fl989l) : 

7* = 8 ■ 10- 10 Tp 6 8 7 cm 3 s- 1 ; 
7 2 e ! = 2.8 • Hr 7 !^' 5 cm 3 s- 1 ; 
A21 = 2.4 ■ 1(T 6 s _1 ; 

AE21 = 1.259 ■ 10" 14 erg (T cxc = 91.2 K, A = 157.74 /jm). 
The cooling function is computed according to formula (23). 

* Sill: we model Sill as a two-level system with the fine 



structure transition (3p)[ 2 P j/ 2 
IHollenbach & McKed fl989h : 



Pi 



/aJ 



Data were taken from 
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(T cxc = 413.6 K, A = 34.8 /an) 
The cooling function is computed following formula ( 1231 ). 



* OI: neutral oxygen is a metastable system formed 
by the (S=l, L=l) triplet and (S=0, L=0,2) doublet, 
(2p)[ 3 P 2 - 3 Pi - 3 P - 1 D 2 - 1 ffo Un order of increasing level , 
with the following ex citations rates (Hollenbach & McKee| [l989l : 
ISantoro & Shuil2 006): 
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The radiative transition probabilities are dOsterbrockl 1 19891: 



Ihe radiative transition p i 
IHollenbach & McKedl 19891) : 
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energy separations are derived from IHollenbach & McKeeld 19891) : 



AP21 = 3.144 ■ 10" 
AP32 = 1.365 ■ 10" 
AP43 = 3.14-10" l; 
AP53 = 3.56 -10 -1: 



14 erg (T cxc = 227.7 K, A = 63.18 /jm); 
14 erg (T cxc = 98.8 K, A = 145.5 fim); 

erg (Tcxc = 2.283 ■ 10 4 K, A = 6300 A); 

erg (Tcxc = 2.578- 10 4 K, A = 5577 A). 



To compute the cooling function, we solve for the five level 
populations and sum over the contributions from each of them. 

* Fell: we adopt a model for a five-level system including the 
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transitions (3d)[ b D g/2 - b D 7/2 - b D 5/2 - 6 B _ 
order of increasing level. For the data see also ISantoro & Shull 
d2006h and references therein: 
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Figure Bl. Scheme of the level models adopted for the different atoms with 
respective line transition data. 
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we assume a fiducial normalization of 10 -5 for missing data on 
e-impact rates. We have checked that the level populations are 
almost insensitive to the adopted values. 
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AE21 = 7.64 ■ 10~ 14 erg 
AE 32 = 5.62 • 10" 14 erg 
AE i3 = 3.87- 10- 14 erg 
AE 5i = 2.27- 10" 14 erg 



(T cxc = 553.58 K, A = 25.99 H 
(T cxc = 407.01 K, A = 35.35 /jm) 
(T cxc = 280.57 K, A = 51.28 /jm) 
(T cxc = 164.60 K, A = 87.41 /jm) 



To compute the cooling function, we solve for the five level 
populations and sum over the contributions from each of them. 



A scheme of the atomic states, with wavelengths of the 
transitions between different levels, is given in Figure IbTI 
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